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turn  and  energy  even  during  grid  restructuring.  Examples  are  taken 
from  simulations  of  shear  flow  and  flow  over  a hydrofoil  in  which  the 
restructuring  algorithms  are  crucial.  Although  the  structure  of  the 
code  is  highly  scalar,  techniques  are  outlined  for  producing  efficient 


code  even  for  the  new  vector  computers. 
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TRANSIENT  FREE-SURFACE  HYDRODYNAMICS 


INTRODUCTION 


I, 


The  hydrodynamics  code  SPLISH  is  designed  for  Lagrangian  simula- 
tions of  transient  free-surface  phenomena.  The  present  version  of  the 
code  was  developed  for  inviscid,  incompressible  flows  in  two  dimensions. 
The  method  uses  a triangular  finite-difference  grid  in  which  triangle 
sides  are  aligned  along  free-surfaces,  interfaces,  boundaries  and  the 

perimeters  of  submerged  bodies.  The  grid  internal  to  these  surfaces 

1 2 

is  left  free  to  reconnect,  adjusting  to  the  time-dependent  flow.  ’ 

In  addition,  vertices  can  be  added  or  subtracted  as  they  accumulate  or 
become  sparse  in  convergent  and  divergent  regions  of  flow.  The  added 
flexibility  gained  through  such  grid  restructuring  permits  the  applica- 
tion of  Lagrangian  techniques  to  large  classes  of  problems  which  were 

formerly  considered  solveable  only  with  the  aid  of  diffusive  Eulerian 
3 4 

rezone  methods.  ’ For  example,  the  simulation  of  shear  flows  and 
flows  about  obstacles  are  possible  with  only  local  changes  in  the  grid. 
This  paper  will  present  the  formulation  and  the  motivation  of  several 
such  grid  restructuring  techniques,  the  algorithms  used  in  implementing 
them  and  examples  of  their  use  in  SPLISH.  Because  the  lack  of  global 
ordering  in  a reconnecting  grid  is  a drawback  to  its  implementation, 
a discussion  of  techniques  to  produce  more  efficient  codes  is  included. 
Examples  will  be  given  of  calculations  performed  on  NRL’s  TI  ASC  pipe- 
line computer. 

Note:  Manuscript  submitted  November  7,  1977. 


THE  CONTROL-VOLUME  APPROACH 


For  a Lagranglan  formulation  of  inviscid,  incompressible  flow, 
Che  basic  equations  are 


and 


VP  - pgy 


(1) 


V . V - 0 


(2) 


together  with  the  conservation  of  mass,  momentum  and  energy.  There  are, 
of  course,  many  possible  ways  to  design  finite-difference  schemes  for 
these  equations.  However,  it  is  in  general  not  possible  to  determine 
which  approach  will  be  successful.^  For  the  case  of  triangular  grids, 
and  in  particular  reconnecting  grids,  there  does  not  exist  a litera- 
ture of  proven  techniques.  Therefore  the  method  chosen  for  SPLISH 
was  the  control-volume  approach,  in  which  the  finite-difference 
equations  are  formulated  to  satisfy  the  conservation  laws  macroscopi- 
cally,  over  a computational  cell.  In  this  way  the  conservation  of 

physical  quantities  is  explicitly  satisfied  by  the  scheme  at  the  outset, 

6 7 8 

and  corrections  for  non-conservation  are  eliminated.  ’ ’ 

The  definition  of  the  control  volume  will  of  course  depend  on 
the  location  at  which  physical  variables  are  applied  on  the  grid.  It 
is  natural  to  specify  positions  and  pressures  at  vertices,  since  the 
Lagranglan  surfaces  coincide  with  surfaces  on  which  pressure  is  defined 
as  a boundary  condition.  In  our  formulation  velocities  and  densities 
are  triangle  centered,  yielding  a staggered  mesh.  Pressure  gradients 
are  therefore  piecewise  linear  within  each  triangle  and  discontinuous 
at  triangle  sides,  as  are  the  triangle  velocities  and  densities. 
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Therefore  all  Che  variables  in  Eq.  (1)  are  triangle  centered  and  it 
is  easy  to  advance  Che  triangle  velocities  either  implicitly  or 
explicitly. 

The  vertex-centered  control  volumes,  CV,  are  used  to  define  the 
new  pressures  through  Eq.  (2),  expressed  as  an  integral  invariant: 

1 V . V dV  » 0 (3) 

0 

CV 

That  is,  the  pressures  at  the  vertices  are  iterated  until  the  resultant 
triangle  velocities  reflect  a divergence-free  condition  for  each  con- 
trol volume.  An  obvious  construction  for  a control  volume  for  this 
application  is  shown  in  Figure  1.  The  vertex-centered  control  volume 
is  defined  by  lines  extending  from  the  triangle  centroids  to  the 
triangle  side  midpoints.  This  permits  a unique,  uniform  and  complete 
tessalation  of  the  entire  computational  region.  The  control  volume 
for  each  vertex  contains  exactly  one-third  of  the  area  of  each  of  the 
adjacent  triangles.  Because  pressures  are  defined  as  boundary  condi- 
tions, the  control  volumes  are  altered  near  boundaries  as  shown  in 
Figure  2.  In  this  way  the  pressures  at  vertices  near  the  boundaries 
enforce  a divergence-free  condition  over  the  additional  area  as  well. 

There  is  another  constraint  implicit  in  Eq.  (1).  Taking  the 
curl  of  both  sides,  we  have 


For  homogeneous  systems  this  implies  that  V x V is  an  invariant  for 
every  control  volume  since  V x 7P  = 0.  With  a properly  defined  7p^ 
the  finite  difference  formulation  also  yields  ^ X 7p  = so  that  the 
vorticity  cannot  be  altered  by  the  pressure  gradients  alone.  However, 
this  invariant  does  not  ensure  the  conservation  of  vorticity  over  a 
complete  timestep  by  itself.  The  velocities  are  triangle-centered,  and 
advancing  the  vertex  positions  means  altering  the  size  and  shape  of  the 
control  volumes  while  leaving  the  velocities  unchanged.  In  other  words 
the  integrated  region  is  changed  but  not  the  integrand.  Therefore  the 
updating  of  vertex  positions  dictates  an  accompanying  change  in  triangle 
velocities  to  keep  the  vorticity  conserved.  This  is  shown  explicitly  in. 
Figure  3.  The  rotation  and  stretching  of  the  triangle  has  necessitated 
a rotation  and  diminishing  of  the  triangle  velocity  such  that  the  V • <il 
line  integral  contributions  within  the  triangle  for  each  of  the  three 
vertex  control  volumes  remains  constant. 

Thus  far  we  have  shown  only  that  a logical  extension  of  the  con- 
trol volume  approach  is  applicable  to  a general  triangular  mesh.  It 
can  lead  to  finite-difference  formulations  which  macroscopically  con- 
serve appropriate  physical  quantities,  regardless  of  how  irregular  the 
mesh  becomes.  The  real  utility  of  this  approach  can  be  seen,  however, 
when  we  allow  the  mesh  the  freedom  to  reconnect. 

RECONNECTIONS 

Despite  the  assurance  that  both  the  curl  and  divergence  can  be 
conserved,  solutions  through  such  Lagrangian  algorithms  can  still 


yield  grossly  erroneous  results.  Figure  4a  illustrates  a sample  cal- 
culation of  shear  flow.  Triangles  below  the  center  of  the  fluid  are 
moving  to  the  left  with  velocity  - U,  and  those  above  are  moving  to  the 
right  with  velocity  U.  The  vertices  lying  directly  on  the  shear  inter- 
face are  stationary,  while  those  above  and  below  move  with  the  triangle 
velocities.  The  shear  layer  is  in  an  unstable  equilibrium  and  should 
persist  until  round-off  errors  accumulate  enough  to  perturb  the  layer. 
However,  even  in  the  absence  of  round-off  in  the  positions,  the  scheme 
will  quickly  fail  since  the  triangles  on  either  side  of  the  boundary 
will  become  stretched.  Very  soon  pressure  gradients  will  be  calculated 
which  involve  vertices  far  removed  from  each  other  as  shown  in  Figure 
4b.  Although  the  area  of  each  triangle  remains  constant  (and  each 
control  volume  divergence-free)  the  convergence  of  the  iterations  would 
slow  and  truncation  errors  build  rapidly  as  the  triangle  sides  lengthened. 
The  very  appearance  of  the  grid  reflects  the  non-physical  situation  in 
which  pressures  far  removed  from  their  co-triangular  points  on  the  in- 
terface directly  influence  their  behavior,  whereas  those  in  the  immediate 
vicinity  have  little  effect. 

Clearly,  allowing  the  mesh  to  reconnect  can  solve  this  dilemma. 

After  reconnection  the  finite-differences  will  again  only  involve 
neighboring  vertices.  The  most  obvious  criterion  for  reconnection  is 
based  on  this  premise.  Any  interior  mesh  line  is  associated  with  a 
triangle  on  either  side.  The  line  can  therefore  be  viewed  as  one  of 
two  possible  diagonals  of  the  quadrilateral  formed  by  these  two 
triangles.  One  reconnection  prescription  is  to  select  the  shorter  of 
the  two  diagonals,  provided  the  resulting  triangles  are  not  too  unequal 
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in  size.  This  last  caveat  can  be  quantified  in  a number  of  ways,  of 
course.  For  example,  this  algorithm  would  remove  the  line  from  vertex 
1 to  vertex  3 and  subtract  a line  from  vertex  2 to  vertex  4 only  in 
Figure  5a.  In  Figure  5b  the  new  line  would  lie  outside  the  quadri- 
lateral and  in  Figure  5c  the  post-reconnection  triangles  would  be 
vastly  different  in  size.  Thus  reconnection  should  not  be  performed. 

Despite  the  simplicity  and  physical  motivation  of  such  an  algorithm, 
it  is  not  obvious  that  it  is  the  preferred  one.  The  expression  for  a 
general  triangular  mesh  Poisson  equation  is 


- - 

X • z = 0 = A (^  • V) 

2 cv  cv 


(5) 


The  vertex  c is  the  central  vertex  as  shown  in  Figure  6a.  The 
is  a sum  over  all  triangles  adjacent  to  c and  the  labelling  for  the 
i + ^th  triangle  is  as  shown  in  the  figure.  is  the  area  of  the 

i + -^th  triangle.  is  the  area  of  the  control  volume  about  vertex 

c and  • V)  is  the  finite  difference  form  of  Eq.  (3).  The  coef- 

CV 

ficient  a of  the  cp  term  is 
c X 


’S) 


kA 


(6) 


and  is  always  negative.  The  coefficient  of  the  cp^  term  has  contribu- 
tions from  just  two  triangles  and  and  is 
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Therefore  the  matrix  represented  by  Eq.  (5)  is  diagonally  dominant  if 
9^  + 9^  fi  l30°  for  each  i.  This  provides  another  uniquely  defined  re- 
connection algorithm,  since  the  sum  of  both  such  pairs  of  angles  in 
the  quadrilateral  is  just  360°.  Whenever  ^ l80°,  the  line  is 

reconnected  to  the  opposite  diagonal.  In  other  words,  we  have  -a  re- 
connection algorithm  which  automatically  ensures  the  diagonal  dominance 
of  the  Poisson  Equation  for  an  irregular  triangular  mesh.  This  al- 
gorithm also  automatically  ensures  that  a diagonal  outside  the  quad- 
rilateral is  never  chosen. 

Whichever  reconnection  algorithm  is  chosen,  during  a reconnection 
the  smallest  physically  definable  cell  is  the  quadrilateral,  and  not 
the  triangles.  It  is  reasonable  then  to  ensure  that  quadrilateral 
properties  are  unchanged  during  a reconnection.  In  other  words,  the 
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quadrilateral  becomes  a control  volume  over  which  certain  variables 
must  be  conserved.  The  reconnection  is  further  complicated  by  its 
alteration  of  th^  vertex  control  volumes  of  each  of  its  vertices. 
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i 

S 

1 
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as  shown  in  ’^igure  7.  Nevertheless,  it  would  be  desirable  that  the 
conservation  laws  enforced  through  vertex  control  volumes  remain  in 
force  during  the  reconnection.  For  example,  to  keep  the  vorticity 

I'  — — 

and  divergence  conserved  the  portions  of  the  integrals  V • V dV  and 
r - - - 

7 X V . dA  within  the  old  and  new  triangles  must  be  the  same 
surface 

before  and  after  the  reconnection.  Since  there  are  four  vertices  with 
two  equations  for  each  and  only  four  triangle  velocity  components  to  be 
changed,  it  would  seem  that  both  integrals  cannot  remain  unaltered. 

In  fact  the  eight  equations  are  not  independent.  There  exists  a unique 
solution  for  every  interior  vertex  which  satisfies  both  constraints. 

It  is  given  by 


where  the  vector  definitions  of  A,  B,  C,  D and  R are  given  in  Figure  8, 
and  are  the  triangle  velocities  before  reconnection  and 
and  Aj,  Ap  are  the  triangle  velocities  and  areas  after  the  reconnec- 
tion. 
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This  solution  includes  a valuable  bonus.  Not  only  is  it  unique, 
but  it  is  reversible.  Re-reconnecting  a line  yields  the  original 


triangle  velocities  as  well  as  the  original  vertex  control  volume  con- 
figurations. This  is  a desirable  feature  since  it  mirrors  yet  another 
implicit  property  of  the  basic  equations,  time  reversibility. 


As  stated  above,  Eq.  (10)  also  preserves  the  quadrilateral  velocity. 


That  is 


Vq  ■ Vr  ■ Vf  * Vs 


where 


The  mass  of  the  quadrilateral  is  just  the  sum  of  the  triangle  masses,  or 


Therefore  strict  mass  conservation  can  also  be  enforced  despite  the 

destruction  and  creation  of  triangles  by  constraining  the  new  densities 

through  Eq.  (13).  By  Eq.  (11)  and  (13)  abd  are  both  separately 

conserved,  and  therefore  so  is  the  quadrilateral  momentum  and  kinetic 

energy.  The  pressures  are  defined  at  vertices  whose  positions  do  not 

change  during  reconnection,  and  the  potential  energy  can  be  altered  only 

V P 

by  the  different  definitions  of  — . Since  we  are  fr“e  to  choose  both 

P 

new  densities,  and  have  only  one  equation,  Eq.  (13),  to  satisfy,  we 
also  control  the  change  in  potential  energy  through  Equation  (4) , which 
provides  the  second  density  constraint.  This  specifies  that  the  amount 


of  vorticity  generated  within  the  quadrilateral  is  the  same  before 
and  after  the  reconnection. 

Reconnections  then  offer  a very  attractive  alternative  to  global 
rezoning.  The  control  volume  approach  leads  to  algorithms  which  conserve 
vorticity  and  divergence  in  the  control  volume  about  each  vertex  despite 
the  reconnections.  The  algorithms  also  conserve  mass,  energy  and  mo- 
mentum on  the  basis  of  triangular  and  quadrilateral  control  volumes  and 
further  exhibit  time  reversibility.  Finally,  the  prescription  for  the 
occurrence  of  reconnections  can  be  chosen  to  preserve  ordering  by 
nearest  neighbors  or  to  preserve  the  diagonal  dominance  of  Poisson's 
equation  over  the  irregular  grid. 

These  routines  have  been  tested  on  the  problem  of  shear  layer  in- 
7 

stability.  Figure  9 illustrates  a simulation  in  which  the  initial 
shear  layer  was  not  perturbed.  The  grid  has  not  changed  although  each 
vertex  in  the  upper  half  of  the  fluid  has  traversed  the  entire  grid  ten 
times.  Figure  10  shows  the  grid  and  a marker  particle  display  for  the 
case  of  an  initial  perturbation  which  has  grown  to  a Kelvin-Helmhoitz 
billow  (top  row),  and  after  the  mature  billow  begins  to  shear.  The 
calculation  agrees  well  with  the  predicted  growth  rates.  Only  the  re- 
connection algorithms  were  used  in  restructuring  the  grid,  and  the  grid 
at  t ■ .329  sec.  exhibits  two  irregularities.  Two  of  the  vertices  have 
become  too  close,  forming  thin  elongated  triangles,  and  a third  vertex 
has  become  enclosed  within  a triangle.  Other  grid  restructuring  is 
clearly  needed.  These  additional  techniques  will  be  discussed  in  the 
following  section. 
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VERTEX  ADDITION  AND  DELETION 


The  fluid  flow  near  a separatrix  is  another  area  in  which  tradi- 
tional numerical  Lagrangian  treatments  fall.  Figure  11  illustrates  a 
sample  grid  for  a submerged  hydrofoil  near  a free  surface.  The  flow 
is  directed  to  the  right  initially  with  velocity  U.  Clearly  as  the  flow 
develops  vertices  will  tend  to  accumulate  at  the  forward  stagnation 
point.  At  the  same  time  the  vertices  on  the  hydrofoil  will  move  with 
the  flow  along  the  hydrofoil  and  accumulate  (by  pairs)  in  its  wake. 

After  a very  short  time  the  gridding  will  deteriorate  to  a denser  grid 
before  and  after  the  hydrofoil  and  a sparse  grid  along  the  hydrofoil 
Itself. 

These  problems  do  not  derive  from  an  ill-chosen  initial  grid  but 
are  inherent  in  Lagrangian  treatments.  By  choosing  the  mid-line  of  the 
hydrofoil  between  an  initial  row  of  vertices,  a situation  develops  in 
which  vertices  of  the  same  triangle  flow  toward  opposite  sides  of  the 
obstacle,  leading  to  a nonsensical  calculation  of  fluid  pressure  gradients 
between  vertices  separated  by  an  obstacle.  Fixing  the  vertices  on  the 
hydrofoil  surface  and  allowing  reconnections  is  neither  physically  de- 
sirable nor  effective.  A tangled  grid  results  from  inverted  triangles 
forming  over  the  hydrofoil.  The  advantage  of  the  present  treatment 
using  a control  volume  formulation  is  that  while  these  problems  are 
difficult,  they  are  still  soluble. 

The  problem  of  representing  flow  near  separatrices  could  be  re- 
solved if  it  were  possible  to  add  and  subtract  vertices  from  the  calcu- 
lation as  needed  without  altering  the  physically  conserved  properties 
of  the  flow.  Fortunately,  such  schemes  are  possible,  and  may  be  derived 


in  the  same  spirit  as  the  reconnection  algorithm  by  using  a control 
volume  approach.  For  example,  the  simplest  way  to  add  a vertex  to  the 
fluid  is  at  the  centroid  of  a triangle.  Lines  are  drawn  from  the  new 
vertex  to  each  of  the  existing  vertices,  leaving  three  new  triangles  in 
place  of  the  old  one.  If  all  triangle  variables  are  chosen  to  be 
identical  with  the  original  variables  of  the  old  triangle,  the  behavior 
of  the  three  triangles  clearly  does  not  alter  the  fluid  motion  since  it 
is  identical  to  that  of  the  initial  triangle.  If  the  vertex  pressure 


is  chosen  as  the  average  of  the  three  vertex  pressures,  we  have  like- 
wise left  unaffected  the  pressure  gradients  throughout  the  area  occupied 
by  the  old  triangle.  In  other  words,  the  three  new  triangles  behave 
exactly  as  the  former  one,  and  are  indistinguishable  from  it  since  the 
vertex  remains  enclosed  within  its  former  boundaries.  However,  we  know 
that  the  reconnection  algorithm  will  eventually  alter  one  of  the  triangle 
sides  if  a vertex  in  the  vicinity  of  the  old  triangle  was  really  needed. 


Otherwise  the  new  vertex  will  continue  to  behave  as  if  it  were  not 
there.  But  the  reconnection  is  also  conservative  as  shown  above,  and 
once  a reconnection  occurs  we  have  successfully  Introduced  a vertex  while 
maintaining  strict  conservation  of  flow  properties. 

The  converse  is  also  true.  If  a vertex  becomes  enclosed  in  a 
triangle,  the  behavior  of  that  triangle  is  not  altered  if  the  vertex  is 
removed  and  the  new  larger  triangle  given  the  velocity 

Vu  = + A^V^  ilk) 


i 

I, 
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where  triangle  4 encompasses  the  three  triangles  1,  2 and  3.  Since  the 
mass  of  the  resultant  triangle  can  be  defined  in  a similar  manner,  the 
momentum  of  the  larger  triangle  has  not  been  altered.  Therefore  for 
both  addition  and  subtraction  of  vertices,  the  larger  triangle  acts  as 
a control  volume  in  exactly  the  same  sense  as  the  quadrilateral  for  the 
reconnections.  What  has  been  lost  is  the  information  about  the  behavior 
of  the  pressure  gradients  and  velocity  gradients  within  the  triangle. 
This  Information  has  been  averaged  out  and  replaced  by  a linear  varia- 
tion across  the  triangle.  All  we  have  suffered  Is  a loss  In  resolution, 
exactly  what  we  set  out  to  do  in  removing  the  vertex. 

This  claim  requires  elaboration.  The  removal  of  a vertex  Implies 
the  alteration  of  four  vertex  control  volumes,  one  of  which  Is  removed, 
and  such  a drastic  change  does  not  seem  consistent  with  a mere  change  in 
resolution.  Figure  12  Illustrates  the  triangles  before  and  after 
vertex  removal.  Before  the  vertex  is  deleted  the  relevant  contribution 
to  the  vorticlty  integrals  about  each  vertex  are 
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After  vertex  4 is  eliminated  its  vorticity  must  be  apportioned  to 
vertices  1 through  3 in  some  manner; 
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Eliminating  S 
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Substituting  - A^V^  + A^V^  + A^V^^  into  Eq.  (17)  yields,  after 

some  algebra. 
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H = A. 


where  since  + A^  + A^  = A^.  Therefore,  conserving 

momentum  exactly  over  the  large  triangle  yields  exactly  the  conserva- 
tion of  vorticity  within  the  affected  vertex  control  volumes.  The 
vorticity  carried  by  the  expunged  vertex  is  apportioned  by  area  to  the 
neighboring  vertices.  If  the  deleted  vertex  lies  close  to  one  of  the 
remaining  vertices,  that  vertex  will  carry  most  of  the  reassigned 
vorticity.  Therefore  the  total  vorticity  is  accounted  for  in  a reason- 
able and  natural  manner,  and  momentum  is  still  conserved.  Since  a 
similar  argument  holds  for  the  divergence  equations,  conservation  of 
flow  variables  is,  demonstrated,  and  a loss  in  resolution  is  the  only 
effect. 

For  the  case  of  an  added  vertex,  the  vertex  control  volume  inte- 
grals are  trivially  left  unchanged,  and  the  added  vertex  initially 
carries  no  vorticity.  Vorticity  can  accumulate  about  the  added  vertex 
only  through  reconnections  with  triangles  having  dissimilar  p.  That 
is,  vorticity  is  generated  only  by  density  gradients,  as  expected. 

Therefore  vertices  can  be  added  and  subtracted  within  triangles 
while  conserving  flow  properties  exactly.  In  both  cases  the  usefulness 
of  this  result  derive  from  the  reconnection  algorithm.  Used  in  tandem 
with  addition  and  deletion  within  triangles  it  provides  a general  al- 
gorthm  for  altering  the  grid  without  disturbing  the  fluid  flow.  It  has 
already  been  shown  that  reconnection  used  after  addition  of  a vertex  at 
the  triangle  centroid  liberates  the  vertex  in  a conservative  manner  and 
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permits  it  to  behave  no  longer  as  the  centroid  of  the  triangle.  The 
process  can  also  be  reversed.  The  reconnection  algorithm  can  be  used  to 
isolate  any  vertex  within  a larger  triangle.  Once  this  is  accomplished 

the  vertex  can  be  averaged  out.  f 

The  utility  of  this  technique  is  not  limited  to  interior  mesh 
points.  Figure  13a  Illustrates  a triangle  at  the  leading  edge  of  a sub- 

I 

merged  body.  The  flow  is  forcing  the  triangle  vertices  on  the  boundary  j 

in  opposite  directions,  and  resolution  of  the  leading  portion  of  the  i 

hydrofoil  is  being  lost.  A point  may  be  added  along  the  body  as  in  13b.  , 

The  result  may  be  viewed  as  the  addition  of  a point  within  a triangle,  ? 

but  one  in  which  only  two  of  the  three  smaller  triangles  survive. 

Figure  14a  illustrates  the  reverse  situation  at  the  rear  of  the  : 1 

i 

hydrofoil.  There  a new  line  is  drawn  to  enclose  the  unwanted  vertex  in  ' 1 

ji 

a triangle.  In  Figure  14b  the  vertex  has  been  removed  from  the  body,  | I 

leaving  the  way  clear  to  add  a vertex  in  the  fluid,  if  needed  to  pre-  | ' 

i 

serve  resolution,  as  in  Figure  14c. 

i J 

The  use  of  the  control  volume  approach  has  therefore  made  possible  ; | 

the  dynamic  addition  and  subtraction  of  vertices  exactly  where  desired  i j 

1 ; 

and  in  a fashion  which  locally  and  globally  conserves  the  properties  of  ' 1 

j i 

the  fluid  flow.  The  combined  use  of  local  resolution  alteration  and  ■ i 

reconnection  algorithms  permits  Lagrangian  calculations  of  extremely  1? 

!] 

complicated  flows. 

EFFICIENCY 


It  is  obvious  that  for  a code  such  as  SPLISH  any  global  ordering 
of  the  arid  would  soon  be  invalidated  by  reconnections  and  by  the 


addition  and  deletion  of  vertices.  Since  standard  fast  Poisson  solvers 
rely  on  such  an  implicit  ordering  among  vertices,  it  is  appropriate  to 
conclude  this  paper  with  some  remarks  on  its  efficiency. 

The  new  generation  of  vector  computers  provides  a good  starting 
point  for  such  a discussion.  The  increased  speed  of  these  computers  is 
attained  in  large  part  by  their  ability  to  perform  quickly  a given  op- 
eration on  large  numbers  of  members  of  an  array,  which  are  preferably 
stored  contiguously  in  memory.  Calculations  on  a vector  computer  are 
therefore  performed  operation  by  operation  for  all  array  members  instead 
of  performing  the  whole  sequence  of  operations  In  turn  for  each  member 
of  an  array.  The  ordering  of  the  vertices  becomes  all  Importnat,  and  a 
highly  disordered  code  such  as  SPLISH  Is  almost  totally  unsuitable  for 
efficient  operation.  However,  despite  such  obvious  problems,  the  en- 
tire SPLISH  code  has  been  optimized  for  the  NHL  TI  ASC  vector  computer. 

Our  first  example  Is  the  reconnection  algorithm  Itself.  The  heart 
of  the  reconnection  algorithm  is  based  on  the  quadrilateral  about  the 
line.  However,  every  reconnection  that  is  performed  redefines  the  quad- 
rilaterals for  each  of  the  four  lines  which  make  up  the  original  quad- 
rilateral. This  situation  is  highly  scalar,  in  that  a single  reconnec- 
tion may  invalidate  the  possible  reconnection  of  four  neighboring  lines. 
Therefore  It  Is  Impossible  to  allow  reconnections  to  proceed  in  parallel, 
and  the  complete  calculation  for  one  must  be  performed  before  the  next 
is  initiated. 

However,  even  in  this  situation  some  increase  in  speed  may  be 
gained  through  efficient  coding.  Clearly  a good  deal  of  the  time  in  the 
reconnection  algorithm  is  spent  in  testing  each  line  for  a possible 


' I 


reconnection.  In  general,  very  few  of  the  lines  reconnect  for  a given 
timestep.  The  flow  is  following  local  streamlines.  Therefore  the  test 
can  be  vectorized  provided  its  output  is  a list  of  lines  which  may  want 
to  be  reconnected.  Each  of  these  (few)  lines  is  then  passed  through  the 
scalar  reconnect  routines,  in  which  they  are  retested  and  the  reconnec- 
tion performed  if  it  is  still  desirable.  An  iteration  through  this 
procedure  may  be  desired,  but  in  most  cases  is  not  necessary  since  most 
reconnections  occur  remote  from  each  other.  In  no  realistic  case  tested 
were  more  than  three  iterations  required  to  reach  the  final  grid. 

The  savings  in  computer  time  of  course  depends  on  the  number  of 
lines  reconnected.  For  roughly  one  percent  of  the  lines  reconnecting 
per  timestep  (an  average  case) , the  vectorized  test  followed  by  scalar 
reconnect  is  ten  times  faster.  The  same  is  true  of  all  the  grid  restruc- 
turing algorithms.  Large  saving  in  time  can  accrue  through  vectorizing 
the  tests  which  must  be  performed  on  every  line  or  triangle.  The  grid 
alterations  must  remain  scalar,  but  are  relatively  few  in  number. 

The  second  example  is  the  solution  of  Poisson's  Equation.  In 
SPLISH  the  pressures  are  adjusted  at  each  vertex  iteratively  to  enforce 
a divergence-free  condition  for  each  vertex  control  volume.  As  shown  by 
Eq.  (5),  coefficients  of  each  term  are  expressed  in  terms  of  the  posi- 
tions of  co-triangular  vertices.  Since  there  is  no  global  ordering, 
such  a calculation  accesses  storage  almost  randomly.  That  is,  the 
code  is  highly  scalar,  and  its  efficiency  on  a vector  machine  is  cor- 
respondingly poor. 

Nevertheless,  it  is  also  possible  to  obtain  vectorized  code  in 
this  situation.  The  solution  is  to  precompute  arrays  which  duplicate 
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Che  position  data  for  each  neighboring  vertex.  These  "alias"  arrays  can 
be  ordered  consecutively  in  core  in  exactly  the  ordering  necessary  for 
efficient  vectorized  code.  Therefore  although  a good  deal  of  extra 
scalar  computation  is  performed,  the  increase  in  speed  obtainable  for 
the  vectorized  code  more  than  compensates  for  the  extra  time.  The  use 
of  "alias"  arrays  has  yielded  decreases  in  computation  time  of  roughly 
a factor  of  three.  Typical  timings  for  the  Poisson  solver  are  now  6.8 
milliseconds  per  iteration  for  120  vertices  or  about  60  microseconds 
per  iteration  per  vertex. 

Although  such  increases  in  speed  are  encouraging,  they  are  not 

the  final  solution.  Large  calculations  will  require  faster  solvers. 

The  most  promising  approach  is  through  the  use  of  direct  solvers,  rather 

than  iterative  ones.  Although  the  matrix  representing  Eq.  (5)  does  not 

exhibit  the  ordering  of  rectangular  meshes,  it  is  nonetheless  sparse. 

Furthermore,  if  the  vertices  are  preordered  by  position,  the  non-zero 

members  will  lie  along  rather  diffuse  bands.  Recently,  there  has  been 

an  increase  in  interest  in  fast  solvers  for  such  banded  matrices  and 

9 10 

several  techniques  look  particularly  encouraging.  ’ The  outlook  is 
very  good.  Not  only  is  a large  class  of  problems  now  amenable  to 
Lagrangian  calculations,  but  also  at  a computational  cost  per  zone 
competitive  with  other  techniques. 
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Figure  1 


Definition  of  a control  volume  about  an  interior  vertex, 
V^.  The  area  of  triangle  j is  apportioned  equally  to  the 
control  volumes  about  V|,  V2  and  V3. 
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Figure  3 

Conservation  of  vorticity  while  advancing  vertex  positions. 
The  triangle  velocity  is  altered  such  that  the  V. dX  line 
integral  contributions  within  the  rotated  and  stretched 
triangle  remain  the  same  for  each  vertex. 
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F igure  4 

An  example  of  stretched  grids  in  Lagrangian  simulations. 
The  initial  grid  has  triangle  velocities  to  the  right  in 
the  upper  half  of  the  fluid  and  to  the  left  in  the  lower 
half  of  the  fluid  (a) . Periodic  boundary  conditions  are 
specified  on  the  sides  of  the  region.  The  grid  very 
quickly  distorts  in  the  vicinity  of  the  shear  layer  (b) . 
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Figure  5 

The  "nearest  neighbor"  reconnection  algorithm.  Lines  are 
reconnected  to  choose  the  shortest  diagonal  of  a quadri- 
lateral (a) . Lines  are  not  reconnected  for  "inverted" 
quadrilaterals  even  if  the  exterior  diagonal  is  shorter 
(b) . Nor  is  the  reconnection  performed  if  the  resultant 
triangles  are  too  dissimilar  in  size  (c) . 
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Figure  6 

The  reconnection  algorithm  to  preserve  diagonal  dominance 
of  the  Poisson  Equation.  The  triangle  ^nd  vertex  labelling 
used  In  the  Poisson  Equation  Is  shown  In  (a) . Figure  (b) 
Indicates  the  angles  9"^  and  9”  used  in  the  reconnection 
test  for  the  line  from  vertex  c to  vertex  1. 


Figure  7 

The  alteration  of  control  volumes  by  a reconnection.  Por' 
tions  of  the  control  volumes  about  vertex  1 and  vertex  2 
are  shown  before  and  after  the  reconnection. 


SHEAR  FLOW  WITH  NO  PERTURBATIONS 


Figure  9 

A test  of  the  reconnection  algorithm  in  SPLISH.  The  grid 
at  various  times  is  presented  for  a simulation  of  an  un- 
perturbed shear  layer.  The  grid  at  t - 0 corresponds  to 
that  in  Figure  4a.  However,  here  the  grid  reconnects  as 
it  stretches,  as  shown  between  t ■ .036  sec  and  t - .040 
sec.  In  this  simulation  8t  ■ .004  sec,  |uj  ■ 5 cm/sec  and 
the  length  of  the  system  is  2 cm.  At  t - 2.00  sec  each 
vertex  in  the  upper  layer  has  passed  each  vertex  in  the 
lower  layer  ten  times.  Any  errors  in  assignment  of 
triangle  velocities  would  have  perturbed  the  unstable 
equilibrium  of  the  layer. 
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Figure  10 

The  grid  and  marker  particles  at  two  times  in  a simulation 
of  a perturbed  shear  layer.  Here  8t  ■ .001  sec,  ju| 

■ 5 cm/sec  and  the  length  of  the  system  is  1 cm.  At  t “ 
.300  sec  the  perturbed  layer  has  grown  to  a mature 
Kelvin-Helmholtz  billow.  Only  reconnections  were  used 
in  restructuring  the  grid.  At  t = .329  sec  the  billow  is 
shearing.  Two  grid  anomalies  are  circled. 


Figure  13 


The  addition  of  a vertex  at  a boundary.  The  vertices  on 
the  boundary  are  moving  along  opposite  sides  of  a sub- 
merged body  (a) , and  resolution  is  lost  for  the  leading 
, edge  of  the  body.  In  (b)  a new  vertex  is  added  on  the 

boundary.  The  old  triangle  is  deleted  and  two  new 
triangles  are  added. 
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Figure  14 


The  deletion  of  a boundary  vertex.  The  flow  is  converging 
at  the  trailing  edge  of  a submerged  body  (a) , resulting  in 
a clustering  of  vertices  and  the  formation  of  long  thin 
triangles.  In  (b)  a vertex  is  removed  by  drawing  a new 
line  to  form  the  enclosing  triangle.  In  (cj  a new  vertex 
is  added  within  the  elongated  triangle  to  preserve 
resolution  of  the  flow  at  the  trailing  edge.  Subsequent 
reconnections  will  remove  the  thin  triangles. 
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